// Fast arctan function for double and quaduple precision
// By Stephen L. Moshier, and conversion to degree math by StellarDX

#include "CSE/Base/CSEBase.h"
#include "CSE/Base/MathFuncs.h"

_CSE_BEGIN

_EXTERN_C

/*
    Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, see
    <https://www.gnu.org/licenses/>.
*/

// Numbers-got nerfed code. (数值受害的代码)

const __Float64 __Arctan128F_table_deg[]
{
    0.0000000000000000000000000000000L,
    7.1250163489017975619533008412068L,
    14.036243467926478582892320159163L,
    20.556045219583464308293612747344L,
    26.565051177077989351572193720453L,
    32.005383208083495560790645750405L,
    36.869897645844021296855612559093L,
    41.185925165709645805088586367179L,
    45.000000000000000000000000000000L,
    48.366460663429801121817059807397L,
    51.340191745909909395994137664827L,
    53.972626614896393257287735737805L,
    56.309932474020213086474505438340L,
    58.392497753751097750225810627096L,
    60.255118703057776265097688282113L,
    61.927513064147042834215359681673L,
    63.434948822922010648427806279547L,
    64.798876354524929311671766058437L,
    66.037511025421816760128170608928L,
    67.166345822082457286036707600070L,
    68.198590513648188229755133913056L,
    69.145541960421653161123687627781L,
    70.016893478100021618547809731703L,
    70.820991974189279890739349667634L,
    71.565051177077989351572193720453L,
    72.255328374943067884132194961601L,
    72.897271030947628039567876821289L,
    73.495638618244980919804984380958L,
    74.054604099077145202342310476739L,
    74.577838681261329296690247935448L,
    75.068582821862447103466463017757L,
    75.529705899934111900861474894266L,
    75.963756532073521417107679840837L,
    76.373005140108460180749059913697L,
    76.759480084812795345292703901156L,
    77.124998440387524356814469831339L,
    77.471192290848489231320126438710L,
    77.799531272619215671464200217594L,
    78.111341960372024856719371525284L,
    78.407824589708929774804754883640L,
    78.690067525979786913525494561660L,
    78.959059819676268005534619751660L,
    79.215702132437400041567946636051L,
    79.460816271371770877270190579774L,
    79.695153531233968054716581161360L,
    79.919402012457678707889722218862L,
    80.134193056915632488796653459530L,
    80.340106921557665756434705963333L,
    80.537677791974382608859929458258L,
    80.727398222799698998435007474346L,
    80.909723079177678217019493295852L,
    81.085073042852133130896449382714L,
    81.253837737444791062798919511021L,
    81.416378519886052409943309282484L,
    81.573030978519327244649302689995L,
    81.724107172924799086782696944762L,
    81.869897645844021296855612559093L,
    82.010673233603103007185159778878L,
    82.146686698021781125766634582435L,
    82.278174199869025289964189251378L,
    82.405356631408555015724136744871L,
    82.528440823407624245370684786921L,
    82.647620640107640261912554879291L,
    82.763077974031995550771223361758L,
    82.874983651098202438046699158793L,
    82.983498255277091249464698997111L,
    83.088772880975318979060980682043L,
    83.190949820386594812540621909569L,
    83.290163192243066861532898701183L,
    83.386539517685243541752932832485L,
    83.480198248343014274598707411820L,
    83.571252251170145199365322143521L,
    83.659808254090090604005862335173L,
    83.745967256083532246832298103540L,
    83.829824904970391135909214452025L,
    83.911471845804825100476820335982L,
    83.990994042505474956721419026891L,
    84.068473075080237482331665150499L,
    84.143986414571043027412823436317L,
    84.217607677635959504969517673186L,
    84.289406862500357487304118651766L,
    84.359450567843177147625729315379L,
    84.427802196036204353752366794796L,
    90.000000000000000000000000000000L
};

const __Float64 __Arctan128F_table_rad[84] =
{
    0.0000000000000000000000000000000000000000E+0L,
    1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125)  */
    2.4497866312686415417208248121127581091414E-1L,
    3.5877067027057222039592006392646049977698E-1L,
    4.6364760900080611621425623146121440202854E-1L,
    5.5859931534356243597150821640166127034645E-1L,
    6.4350110879328438680280922871732263804151E-1L,
    7.1882999962162450541701415152590465395142E-1L,
    7.8539816339744830961566084581987572104929E-1L,
    8.4415398611317100251784414827164750652594E-1L,
    8.9605538457134395617480071802993782702458E-1L,
    9.4200004037946366473793717053459358607166E-1L,
    9.8279372324732906798571061101466601449688E-1L,
    1.0191413442663497346383429170230636487744E+0L,
    1.0516502125483736674598673120862998296302E+0L,
    1.0808390005411683108871567292171998202703E+0L,
    1.1071487177940905030170654601785370400700E+0L,
    1.1309537439791604464709335155363278047493E+0L,
    1.1525719972156675180401498626127513797495E+0L,
    1.1722738811284763866005949441337046149712E+0L,
    1.1902899496825317329277337748293183376012E+0L,
    1.2068173702852525303955115800565576303133E+0L,
    1.2220253232109896370417417439225704908830E+0L,
    1.2360594894780819419094519711090786987027E+0L,
    1.2490457723982544258299170772810901230778E+0L,
    1.2610933822524404193139408812473357720101E+0L,
    1.2722973952087173412961937498224804940684E+0L,
    1.2827408797442707473628852511364955306249E+0L,
    1.2924966677897852679030914214070816845853E+0L,
    1.3016288340091961438047858503666855921414E+0L,
    1.3101939350475556342564376891719053122733E+0L,
    1.3182420510168370498593302023271362531155E+0L,
    1.3258176636680324650592392104284756311844E+0L,
    1.3329603993374458675538498697331558093700E+0L,
    1.3397056595989995393283037525895557411039E+0L,
    1.3460851583802539310489409282517796256512E+0L,
    1.3521273809209546571891479413898128509842E+0L,
    1.3578579772154994751124898859640585287459E+0L,
    1.3633001003596939542892985278250991189943E+0L,
    1.3684746984165928776366381936948529556191E+0L,
    1.3734007669450158608612719264449611486510E+0L,
    1.3780955681325110444536609641291551522494E+0L,
    1.3825748214901258580599674177685685125566E+0L,
    1.3868528702577214543289381097042486034883E+0L,
    1.3909428270024183486427686943836432060856E+0L,
    1.3948567013423687823948122092044222644895E+0L,
    1.3986055122719575950126700816114282335732E+0L,
    1.4021993871854670105330304794336492676944E+0L,
    1.4056476493802697809521934019958079881002E+0L,
    1.4089588955564736949699075250792569287156E+0L,
    1.4121410646084952153676136718584891599630E+0L,
    1.4152014988178669079462550975833894394929E+0L,
    1.4181469983996314594038603039700989523716E+0L,
    1.4209838702219992566633046424614466661176E+0L,
    1.4237179714064941189018190466107297503086E+0L,
    1.4263547484202526397918060597281265695725E+0L,
    1.4288992721907326964184700745371983590908E+0L,
    1.4313562697035588982240194668401779312122E+0L,
    1.4337301524847089866404719096698873648610E+0L,
    1.4360250423171655234964275337155008780675E+0L,
    1.4382447944982225979614042479354815855386E+0L,
    1.4403930189057632173997301031392126865694E+0L,
    1.4424730991091018200252920599377292525125E+0L,
    1.4444882097316563655148453598508037025938E+0L,
    1.4464413322481351841999668424758804165254E+0L,
    1.4483352693775551917970437843145232637695E+0L,
    1.4501726582147939000905940595923466567576E+0L,
    1.4519559822271314199339700039142990228105E+0L,
    1.4536875822280323362423034480994649820285E+0L,
    1.4553696664279718992423082296859928222270E+0L,
    1.4570043196511885530074841089245667532358E+0L,
    1.4585935117976422128825857356750737658039E+0L,
    1.4601391056210009726721818194296893361233E+0L,
    1.4616428638860188872060496086383008594310E+0L,
    1.4631064559620759326975975316301202111560E+0L,
    1.4645314639038178118428450961503371619177E+0L,
    1.4659193880646627234129855241049975398470E+0L,
    1.4672716522843522691530527207287398276197E+0L,
    1.4685896086876430842559640450619880951144E+0L,
    1.4698745421276027686510391411132998919794E+0L,
    1.4711276743037345918528755717617308518553E+0L,
    1.4723501675822635384916444186631899205983E+0L,
    1.4735431285433308455179928682541563973416E+0L, /* arctan(10.25) */
    1.5707963267948966192313216916397514420986E+0L  /* pi/2 */
};

/*
 *	Inverse circular tangent for quaduple precision based on degrees
 *      (arctangent)
 *
 * DESCRIPTION:
 *
 * Returns degree angle between -90 and +90 whose tangent is x.
 *
 * The function uses a rational approximation of the form
 * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
 *
 * The argument is reduced using the identity
 *    arctan x - arctan u  =  arctan ((x-u)/(1 + ux))
 * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
 * Use of the table improves the execution speed of the routine.
 *
 * ACCURACY:
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -19, 19       4e5       1.7e-34     5.4e-35
 *
 * WARNING:
 *
 * This program uses integer operations on bit fields of floating-point
 * numbers.  It does not work with data structures other than the
 * structure assumed.
 */

Angle __cdecl __IEEE754_ATAN128F_C64F(__Float64 x)
{
    int32_t k, sign, lx;
    __Float64 t, u, p, q;
    __Float64 xhi;

    xhi = x;
    k = xhi.parts.msw;
    lx = xhi.parts.lsw;
    sign = k & 0x80000000;

    #ifndef TRIGONOMETRY_USE_RADIANS
    const __Float64* __ArctanF128_table = __Arctan128F_table_deg;
    #else
    const __Float64* __ArctanF128_table = __Arctan128F_table_rad;
    #endif

    /* Check for IEEE special cases.  */
    k &= 0x7fffffff;
    if (k >= 0x7ff00000)
    {
        /* NaN. */
        if (((k - 0x7ff00000) | lx) != 0) { return __Float64::FromBytes(BIG_NAN_DOUBLE).x; }

        /* Infinity. */
        if (sign) { return -__ArctanF128_table[83].x; }
        else { return __ArctanF128_table[83].x; }
    }

    if (k <= 0x3c800000) /* |x| <= 2**-55.  */
    {
        /* Raise inexact.  */
        if (1e300L + x > 0.0) { return x.x; }
    }

    if (k >= 0x46c00000) /* |x| >= 2**109.  */
    {
        /* Saturate result to {-,+}90.  */
        if (sign) { return -__ArctanF128_table[83].x; }
        else { return __ArctanF128_table[83].x; }
    }

    if (sign) { x = -x; }

    if (k >= 0x40248000) /* 10.25 */
    {
        k = 83;
        t = -1.0 / x;
    }
    else
    {
        /* Index of nearest table element.
       Roundoff to integer is asymmetrical to avoid cancellation when t < 0
           (cf. fdlibm). */
        k = int32_t(8.0 * x + 0.25);
        u = 0.125 * k;
        /* Small arctan argument.  */
        t = (x - u) / (1.0 + x * u);
    }

    /* Arctan of small argument t.  */
    static __Float64
        p0 = -4.283708356338736809269381409828726405572E+1L,
        p1 = -8.636132499244548540964557273544599863825E+1L,
        p2 = -5.713554848244551350855604111031839613216E+1L,
        p3 = -1.371405711877433266573835355036413750118E+1L,
        p4 = -8.638214309119210906997318946650189640184E-1L,
        q0 = +1.285112506901621042780814422948906537959E+2L,
        q1 = +3.361907253914337187957855834229672347089E+2L,
        q2 = +3.180448303864130128268191635189365331680E+2L,
        q3 = +1.307244136980865800160844625025280344686E+2L,
        q4 = +2.173623741810414221251136181221172551416E+1L;

    u = t * t;
    p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
    q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
    u = t * u * p / q + t;

    /* arctan x = arctan u  +  arctan t */
    auto __ConvertUnit = [](float64 x) {return Angle::FromRadians(x);};
    u = __ArctanF128_table[k].x + __ConvertUnit(u);
    if (sign) { return (-u.x); }
    else { return u.x; }
}

_END_EXTERN_C

_CSE_END
